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1.  Introduction 

Let  X  =  (Xt:t  >0}  be  a  regenerative  process  which  we  wish  to 
simulate.  Under  mild  regularity  conditions  the  distribution  of  Xt  con¬ 
verges  to  the  distribution  of  some  limiting  random  variable  (or  vector)  X. 
This  type  of  convergence  is  known  as  weak  convergence  and  written 
Xt  =^X}  as  t  t  oo.  Simulators  speak  of  X  as  the  "steady- state"  con¬ 
figuration  of  the  system  and  are  often  interested  in  estimating  the 
constant  r  =  E(f(X)),  where  f  is  a  given  real-valued  function  defined 
on  the  state- space  of  the  process  X.  The  regenerative  method  of  estima¬ 
tion  provides  a  means  of  constructing  point  and  interval  estimates  for  r; 
see  IGLEHART  (1977)  for  an  expository  summary  of  this  method. 

The  problem  we  consider  in  this  paper  does  not  involve  estimation 
of  r,  but  rather  the  estimation  of  extreme  values  of  the  regenerative 
process  X.  Suppose,  for  the  sake  of  discussion,  we  are  simulating  a 

stable  GI/G/l  queue  in  order  to  estimate  the  maximum  waiting  time  among 

* 

the  first  n+1  customers;  call  this  random  variable  As  n  grows, 

so  will  W  .  However,  W  does  not  converge  to  a  finite  limit,  but  rather 
n  ’  n  ' 


diverges  to  +«.  We  will  be  interested  in  estimating  the  distribution 
function  of  Wr  for  finite,  but  large  n.  By  the  same  token  we  might 
wish  to  estimate  the  distribution  of  the  maximum  queue  length  during  the 

I 

interval  [0,t].  While  this  problem  of  estimating  extreme  values  would  seem 
to  be  of  great  practical  importance  to  simulators,  we  know  of  no  papers  in 
the  simulation  literature  which  offer  any  guidance  on  the  subject.  This 
paper  will  attempt  to  partially  fill  the  gap. 

We  begin  in  Section  2  by  summarizing  a  series  of  probabilistic 
results  in  extreme  value  theory  which  will  provide  the  theoretical  basis  for 
the  methods  we  propose.  Section  3  discusses  several  methods  for  estimating 
extreme  values  for  the  general  regenerative  simulation.  In  Sections  lv  and  5 
we  treat  the  special  cases  of  the  GI/G/1  queue  and  birth-death  processes 
respectively.  Theoretical  results  are  available  for  these  two  classes 
of  regenerative  processes  that  are  useful  in  assessing  the  accuracy  of 
the  estimation  methods  proposed.  Section  6  contains  the  numerical  results 
for  simulations  of  the  M/M/ 1  queue  carried  out  to  illustrate  the  estima¬ 
tion  methods  proposed. 


2.  Probabilistic  Background 

Let  {F^rn  >  1)  be  a  sequence  of  distribution  functions  (d.f. 's) 

on  the  real  line, IR  =  (-»,  +m) .  This  sequence  converges  weakly  to  a 

d.f.  G  if  lim  F  (x)  =  G(x)  for  all  x  £]R  which  are  continuity 
n  — *  oo 

points  of  G.  We  write  F  =>  G  to  denote  this  type  of  convergence.  If 

n 

(resp.  X)  is  a  random  variable  (r.v.)  with  d.f.  F^  (reap.  G),  we 
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also  write  X^  ^  X  to  denote  this  weak  convergence.  Sometimes  it  is 

convenient  to  write  X^  =>  G  to  denote  the  same  thing.  The  material 

presented  in  this  section  can  be  £ound  for  the  most  part  in  deHAAN  (1970), 

currently  the  best  comprehensive  treatment  of  the  subject. 

Now  let  (Xn:n  >1]  be  a  sequence  of  independent,  identically 

distributed  (i.i.d.)  r.v.'s  and  denote  the  maximum  of  the  first  n  r.v.'s 

by  M  s  max{X,  :1  <  j  <  n).  If  each  of  the  X  's  has  d.f.  F.  then  M 
n  j  -  J  -  j  >  n 

will  have  d. f.  Fn.  We  shall  say  that  F  belongs  to  the  domain  of 

attraction  of  the  nondegenerate  d.f.  G,  and  write  F  £  £l(G),  if  we  can 

choose  two  sequences  of  constants  {a  :n  >  1)  and  {b  :n  >  1)  with 

n  —  n  — 

a  >0  such  that 

n  • 

(2.1)  Fn(anx  +  bn)  -*G(x) 

as  n  -» oo  for  all  x  £  1R  for  which  G  is  continuous.  Equivalently, 

F  £  &(G)  if  (**n“bn)/an  ^  G  as  n  -» oo.  Thus  for  large  n  we  would 
approximate  P(Mn  <  x)  by  G( ( x-bn)/an) .  If  a  r.v.  X  has  d.f. 

F  £  JS(G),  we  also  write  X  £  fl(G). 

A  famous  result  in  extreme  value  theory  states  that  the  only  d.f. 's 
G  which  can  arise  in  (2.1)  are  of  one  of  the  following  three  types: 

SO  ,  x  <  0 

exp( -x'0)  ,  x  >  0 
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vCt 


exp{-(-x)  }  ,  x  <  0 


(2.3) 


v«)  - 


x  >  0 


(2.4) 


A(x)  =  exp(-e'x)  ,  x  £]R  , 


where  in  (2.2)  and  (2.3)  a  is  a  positive  constant.  Recall  that  two 
d.f.  's  and  Gg  are  said  to  be  of  the  same  type  if  there  exist 

two  constants  a  and  b,  a  >  0,  such  that  G^x)  =  Gg(ax  +  b)  for  all 
x  €1R.  Thus  aside  from  translations  and  scaling  by  a  positive  constant 
the  three  d.f.  's  given  in  (2.2)  -  (2.1+)  are  the  only  ones  that  --in  appear 
in  (2.1).  This  result  on  the  three  types  of  limit  d.f. 's  is  usually 
attributed  to  GNEDENKO  (1943),  however  it  was  first  formulated  in  this 
way  by  FISHER  and  TIPPETT  ( I928) . 

The  next  logical  result  to  seek  is  necessary  and  sufficient  conditions 
for  F  €  fl(G),  where  G  of  necessity  is  one  of  the  three  d.f. 's  given  in 
(2.2)  -  (2.4).  Furthermore,  if  F  E  ■e(G)  we  need  a  method  for  selecting 
the  two  sequences  (a^m  >  1)  and  (b^in  >  lj.  To  this  end  we  first 
define  the  right  endpoint,  x^  <  +«,  of  the  d.f.  F  as 


xQ  =  sup{x:F(x)  <  1)  . 


A  d.f.  F  £  &(  4>a) 
(2.5) 


if  and  only  if  for  all  x  >  0 


lim 

t  — » 00 


1  -  F(tx) 

r-  w 


-a 
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rr 


If  F  €  fi(®a),  c^en  we  can  take 

(2.6)  an  *  lnf(x:l  -  F(x)  <  l/°) 


and  bR  ss  0.  A  d.f.  F  £  if  and  only  if 

x  >  0 


(2.7) 


lim 
t  — *  00 


1  -  F[xQ  -  (tx)'1] 
1  -  F(x0-t_1) 


-a 

X 


If  F  €  j0(¥a)j  Chen  we  can  take  bR  =  xQ  and 


a^  =  Xq  -  inf(x:l  -  F(x)  <  I/n) 


The  final  case,  F  €  J0( A) ,  is  the  most  important 
applications.  A  d. f .  F  €  $(A)  if  and  only  if 

(c.a; 

where  for 


If  F  €  fi(A),  then  we  can  take 


lim  1  -  F(  t  +  xf(  t) )  _  -x 
t  t  xn  1  -  F(t)  =  6 


t  <  x„ 


/  (1  -  F(s))ds 

f<‘>  -  —  r.-fra — 


xQ  <  oo  and  for  all 


one  for  our  simulation 


for  all  x  £]R  , 
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inf(x:l  -  F(x)  <  1/n) 


b  = 
n 


and 


a 

n 


/ 

b 

n 


(1  -  r(t)]dt 


Alternative  expressions  are  available  for  a  and  b  .  We  can  use  a' 

n  n  n 

and  b'  provided  a  /a’  1  and  (b  -b')/a  -+  0.  Let  Q  (p)  denote 

n  n  n  n.  n  n  n 

the  p-quantile  of  the  d.f.  Fn:  for  0  <  p  <  1, 


Qn(p)  =  inf{x:Fn(x)  >  p)  . 
Then  if  F  £  £(A),  we  can  alternatively  select 


and 

(2.1°)  an  "  Ve’e'l>  *  Ve‘l)  • 

Furthermore,  if  F.  £  £(A)  and  xQ  =  +«,  Mn/bn  =*  1  as  n  -*  oo.  Many 
of  the  classical  d.f. 's  such  as  the  exponential,  gamma,  normal,  lognormal, 
and  logistic  belong  to  /£{A) . 

Suppose  F  has  x^  =  +»  and  possesses  an  exponential  tail: 

(2.11)  1  -  F(x)  ~  b  exp(-ax)  ,  as  x  -* «  , 
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where  a  and  b  are  two  positive  constants.  Then  it  is  easy  to  check 

that  (2.8)  holds  and  F  €  S(A).  Using  the  expressions  (2.9)  and  (2.10) 

it  can  be  shown  that  b  and  a  can  be  selected  as  follows: 

n  n 

(2.11a)  bn  =  a_1  ln(tlb)  ' 

and 

-1 

(2.11b)  an  =  a 

An  interesting  (and  practical)  situation  arises  if  F  is  a  discrete 
d.f.  as,  for  example,  the  geometric  d. f .  F(x)  =  1  -  exp(-[x]),  x  >  0, 
where  [x]  is  the  integer  part  of  x.  In  this  case  neither  (2.5^  nor 
(2.8)  hold,  and  since  xQ  =  +»,  F  does  not  belong  to  the  domain  of 
attraction  of  any  of  the  three  types  (2.2)  -  (2.4).  However,  a  result 
has  been  salvaged  by  ANDERSON  ( I97O) .  Let  a  be  the  class  of  all  d.f. 's 
whose  support  consists  of  all  sufficiently  large  positive  integers.  Then 
for  F  €  Q, 


(2.12) 

n  1 

lim  sup  F  (cf  x  +  b  )  <  exp(-e  ) 

n  00 

and 

(2.13) 

lim  inf  Fn(cflx  +  bn)  >  exp(-e‘^X' 

n  a> 


for  some  a  >  0,  all  x,  and  some  sequence  (bn:n  >  1)  if  and  only  if 


(2.14) 


lim 
n  w 


1  -  F(n) 

1  -  F(  n+l)  = 


a 

e 
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When  this  condition  holds,  the  constants  b  can  be  selected  as  follows. 

*  n 

For  F  £  Q,  and  each  positive  integer  n  let  h(n)  =  -log(l-F(n))  and 

define  h  to  be  the  extension  of  h  obtained  by  linear  interpolation 
i  c 

for  x  >  1.  Then  define  for  x  >  1 

Fc(x)  =  1  -  exp(-hc(x))  . 

Clearly  Fc  is  a  continuous  d. f.  and  for  sufficiently  large  x  is 
strictly  increasing  since  F  £  A.  For  x  <  1  the  F£  can  be  defined 
arbitrarily  just  so  long  as  it  is  a  d. f .  In  terms  of  Fc  we  can  define 
bn  for  large  n  as  the  unique  root  of 

1  .  Fe(bn)  .1/0. 

If  F  £  G  and  1  -  F(n)  ~  b  exp(-an)  as  n  -♦  oo  (a^b  >  0),  then 
for  =  a  ^  ln(nb)  it  can  easily  be  shown  using  the  method  followed 
by  HEYDE  ( I97O)  that  for  integer  l 

(2.15)  lim  [?{Mn  -  [bj  <  i}  -  exp(-e'a^-dn))l  =  0  f 

n  -»  00  L  4 

where  d  =  b  -  [b  ].  Thus  for  n  large  we  would  (Ignoring  a  possible 
n  n  n 

continuity  correction)  use  the  approximation 


8 


P{Mn  <  i  +  [bn]  }  a  exp(-e-a(S-dn>)  , 
or 

(2.16)  P(Mn  <  £)  =  exp(-e‘a(A’bn>)  . 

Suppose  now  that  we  also  have  defined  on  the  probability  triple 
( 3,P)  that  supports  the  i.i.d.  sequence  (X^:n  >  1)  a  renewal  process 
{£(t):t  >  0)  with  mean  time  between  renewals  m  (0  <  m  <  «).  Then 
the  weak  law  of  large  numbers  for  renewal  processes  states  that 
:'t)/  t  m  *  as  n  -» oo.  Next  set 

=  max{2$j  :  1  <  j  <  £(  t) )  .  • 

The  following  useful  result  for  this  situation  was  obtained  by  BERMAN 
(I962).  If  (Mn-bn)/an  one  of  the  three  extreme  value  d.f. 's 

(2.2)  -  (2.4),  then  as  t  -*  00 

(2.17)  (M*-b[t])/a[t]  =>Gl/m  . 

This  result  provides  a  useful  tool  for  extreme  values  of  regenerative 
processes.  To  be  explicit  suppose  X  =  (Xt:t  >  0)  is  a  regenerative 
process  defined  on  (fi,3,P)  and  T^ ,  j  >1,  is  the  time  of  the jth 
regeneration  point  of  X  with  =  0.  Then  the  renewal  process 
(£(t):t  >  0)  which  counts  the  number  of  regeneration  points  in  (0,t] 
is  defined  by 
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£(  t)  =  max(j  :T^  <  t} 
with  4(0)  =  0.  For  j  >  1,  let 

I 

=  sup(Xt:Tj_1  <  t  <  Tj)  . 

Since  X  is  re  generative,  the  sequence  of  maxima^  (M.J:j  >  1),  will  be 

i.i.d.  Then  if  L.  =  sup(X  :0  <  s  <  t]  clearly 
t  s  —  — 

(2.18)  max{Mj : 1  <  j  <  £( t) }  <  Lt  <  max(M*:l  <  j  <  4(  t)  +  1)  . 

Combining  the  inequalities  of  (2.18)  with  the  limit  theorem  of  (2.17) 
enables  us  to  show  that 

(2.19)  <VV)>/*[ t)  g1/"  > 

where  m  =  EtT^},  provided  M*  £  )B(G).  Of  course  if  £  Ci,  then  the 
weaker  results  of  Anderson  or  Heyde  are  all  that  can  be  expected. 

We  conclude  this  section  by  summarizing  the  problems  confronting 
us  for  a  regenerative  processes  with  continuous  state  space.  If  M*  £  fl(G), 
then  we  can  use  (2.19)  to  obtain  the  asymptotic  (for  large  t)  approxima¬ 
tion 

(2.20)  P(Lt  <  x)  =  G^"^  a^1)  * 
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Alternatively,  we  can  also  show  that 


(MUt)  '  bJl(t))/aUt)  **  G 

which  when  combined  with  (2.18)  yields  the  asymptotic  (for  large  t) 
approximation 

(2.21)  P{L  <  x}  ^  g( — -  .  . 

c  ”  “  l  I 

This  expression  does  not  require  an  estimate  for  m  as  is  the  case 
with  (2.20). 


If  the  simulation  is  run  for  n  cycles,  then  (2.20)  or  (2.21)  should  be 
replaced  by  • 


(2.22) 


1 1  max  Mt  <  x  l  £  G  ( 
»1<J <n  V 


x-b 


For  (2.20),  (2.21),  or  (2.22)  to  be  useful,  we  must  estimate  an>  and  bQ. 

To  use  (2.20)  we  must  also  estimate  m.  The  expected  cycle  length,  m,  can 
of  course  be  estimated  by  the  sample  mean  of  the  cycle  lengths  observed.  Sev¬ 
eral  methods  for  estimating  a  and  b  will  be  discussed  in  Section  3. 

n  n 

Finally,  we  must  assess  whether  £  <f(G)  for  one  of  the  three  d.f.'s 
(G's)  given  in  (2.2)  -  (2.4).  For  many  simulations  in  which  extreme  values 
are  being  estimated,  the  limit  d.f.'s  G  will  be  either  \  or  $  ,  since  the 
maxima  arising  are  unbounded  (x^  *  +°>) .  Our  experience  with  specific  examples 
indicates  that  if  the  regenerative  process  is  stable  (converges  to  a  non¬ 
degenerate  limit),  then  G  ■  A;  while  if  the  process  is  "null-recurrent" 

(m  =  E(T^)  =  +»),  then  G  =  0^.  However,  in  this  case  G^m(x)  =  1 
for  all  x  >0  which  indicates  that  a  different  normalization  must  be 


used  to  obtain  a  nondegenerate  limit.  In  any  case,  we  note  that  if 

X  £  with  constants  a^  >  0  and  =  0,  then  In  X  £  fi(A)  with 

constants  a  =  <2  and'  b  =  In  a1  .  We  note  in  passing  that  the  extreme 
n  n  n  r  a 

value  behavior  of  some  functions  of  a  regenerative  process  can  be  handled 

I 

in  the  same  way.  If  the  state  space  of  the  regenerative  process  is  discrete, 
then  we  shall  only  consider  the  situation  in  which  the  d.f.  of  M*  £  Q  and 

P(M*  >  n}  ~  b  exp(-an)  as  n  -*■  «  for  some  a,b  >  0.  In  this  case  we  can 
approximate  the  d.f.  of  Lt  or  max{M|:l  <  j  <  nl  by  using  (2.12)  and  (2.13) 
or  (2.15)  and  (2.16). 


3 .  Statistical  Estimation  Problem 

We  now  present  some  procedures  which  can  be  used  to  estimate  the  d.f. 
of  extreme  values  occurring  in  a  regenerative  stochastic  process  belonging 
to  a  (known)  domain  of  attraction.  (The  methods  will  be  Illustrated  for 
the  case  of  maximum  values,  but  analogous  procedures  could  be  used  for 
minimum  values.)  The  basic  idea  underlying  all  the  procedures  is  to 
simulate  a  process  in  order  to  form  an  empirical  d.f.  of  say  the  maximum  of 
the  process  over  a  given  period  of  time  or  number  of  cycles,  and  to  regress 
this  against  the  functional  form  of  the  appropriate  extreme  value 
distribution.  The  methods  differ  as  to  how  the  empirical  d.f.  is  formed 
and  how  the  regression  is  done.  The  applicability  of  the  methods  depends 
on  the  process  lying  in  the  domain  of  attraction  of  a  nondegenerate  d.f. 

For  a  process  taking  on  positive  values,  but  having  right  endpoint  x^  <  «■>, 
convergence  to  a  nondegenerate  extreme  value  distribution  would  not  hold. 
This  would  be  the  case  for  instance  for  the  queue  length  process  in  a 
system  with  finite  waiting  room. 
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3.1.  The  Continuous  State  Space  Case 

Suppose  we  are  Interested  in  the  maximum  of  a  continuous  state  space 
regenerative  process  over  n  cycles.  Choose  an  integer  k  >  n,  (the 
choice  of  k  will  be  discussed  later)  and  simulate  v  cycles  of  the  pro¬ 
cess,  yielding  individual  cycle  maxima  y.  .>  y_  .>...>  y  .  A  sample 

JLfiv  4 

of  k  cycles  contains  (^)  subsets  of  n  cycles.  We  then  use  the 
simulation  results  to  find  the  maximum  of  the  process  over  each  subset  of 
n  cycles.  The  only  y^  k's  which  can  be  maxima  for  some  subset  of  n 
cycles  are  y^  ^ , . . . , ^ ( ^  As  explained  below,  we  are  only  interested 
in  yj^  k,...,yi  k  where  the  predetermined  number  iQ  is  usually 
considerably  less  than  k-n+1.  Thus  in  performing  the  simulation,  we  need 
'•rsly  keep  track  of  y  .  ,  ...,y.  .  ,  and  the  continuous  state  space  assump- 

1  pK.  1q»K. 

tion  can  be  relaxed  to  the  distinctness  of  y,  , . . . ,y  .  (For  instance 

1  f  k  1q  ,  K. 

the  customer  waiting  time  process  of  a  queueing  system  may  have  an  atom  at 

0;  however  this  usually  need  not  concern  us.)  Now  observe  that  y  is 

I 

the  maximum  value  in  (n_^)  the  (^)  subsets  of  n  cycles, 

i  “  1 . k-n+1.  Thus  y^  ^  is  the  maximum  value  in  a  fraction  n/k 

of  the  (n)  cycles,  y ^  k  is  the  maximum  value  in  a  fraction 
n(k-n)/(k(k-l))  of  the  (^)  cycles,  etc.  Let  En(x)  (suppressing  k) 
be  the  empirical  d.f.  of  the  maximum  of  the  process  over  n  cycles.  Then 


En(?l,k> 

-  1, 

V*2,k> 

_  1  /,  k~U 

1  —  n/k  “  k  t 

En^3,k> 

_  k-n  _  n(k-n)  _  k-n 

k-n-1 

k  k(k-l)  k 

’  k-1  » 

'  \-l« 

for  i  ■ 

.  ,ir 


En^i  k^  can  Ee  Stained  ky  a  simple  recursive  computation. 


Thus 


It  is  desirable  to  have  the  values  E^y^  k)  diminish  slowly  enough 
from  1  to  0  so  as  to  provide  good  upper  quantile  information  for  the  re¬ 
gression  procedure.  Thus  we  want  a  rather  high  k/n  ratio  (say,  at  least 

I 

10).  But  n  should  be  large  enough  so  that  the  extreme  value  distribution 
provides  a  good  approximation  to  the  maximum  of  the  process  over  n 
cycles.  So  is  very  small  for  many  of  the  latter  values  of  the 

k-n+1  possible  i's.  Using  too  many  y  's  in  the  regression  would  tend 
to  drown  out  the  effect  of  the  first  few,  and  could  be  a  source  of  noise 
when  E  (y.  ,  )  is  very  close  to  0.  Thus  we  only  consider  y  for 

II  1  )  iC  1  j  K 

i  <  i..,  where  in  -  max{i  :  E  (y,  .  )  >  c)  for  e  such  as  .001. 

—  O’  0  nwi,k  — 

Consider  the  nonlinear  regression  problem  of  selecting  a  -  and  b 

n  n 

so  as  to  obtain  the  least-squares  fit  of  A((Xj^  -  bn)/aQ)  to  En^xi^»  w^en 

A  is  the  appropriate  domain  of  attraction,  for  data  points  xt  **  y2  k . 

y  .  For  short  we  shall  denote  this  problem  by 
*0, 


E  (x. )  -  exp(-exp(-(x  -  b  )/a  )). 
n  1  inn 


When  is  the  appropriate  domain  of  attraction,  the  problem  becomes 


.-a. 


En(Xi)  “  expt-^/a^)  ). 


(Note  that  the  data  form  a  dependent  sample.)  The  point  y  is  not  used 

in  the  regression  since  En(y^  *  1*  which  effectively  makes  y^  ^  the 
right  endpoint  x^,  even  though  A  has  a  right  endpoint  of  +  •  .  This 
manifests  itself  in  the  regression  by  forcing  b^  to  •  or  a^  to  0. 
Refer  to  the  above  procedure  (for  a  fixed  k,n,  and  Iq(O)  as  the  basic 


I 


nonlinear  regression  procedure.  In  the  case  of  domain  A,  this  procedure 

A  A 

provides  estimates  a  ,  b  of  the  norming  constants  for  A.  An  alterna- 

n  n 

tive  procedure,  referred  to  as  the  basic  linear  regression  procedure,  is  to 
(take  a  double  -In  transformation  and)  perform  the  linear  regression 

A 

-ln(-ln(E  (x.)))  *  c  x.  +  d  ,  with  x.  as  above,  to  obtain  estimates  c 
'  n  i  n  1  n  1  n 

A 

and  d  .  This  fit  can  be  used  directly,  or  equivalently  note  that  we  can 
n 

form  estimates  a  ■  c  £  •  -c  ^3  of  the  norming  constants  for  A. 

n  n  n  n  n 

Recall  that  if  X  e  £'($  )  with  parameters  a’  >  0  (and  b’  ■  0) 

a  n  n 

then  In  X  £  with  norming  constants  a  ■  a  ^  and  b  •  In  a*  . 

n  n  n 

Thus  performing  the  nonlinear  regression  procedure  on  X,  based  on 

X  e  A®a),  would  lead  to  E^x^  •  exp^x^/a^j*)  in  order  to  estimate 

a',  a.  Performing  the  nonlinear  regression  procedure  on  In  X,  based  on 
Q 

In  X  eAA),  would  lead  to  EQ(ln  -  exp(-exp(-(ln  Xj-b^/a^))  in  order 

to  estimate  a  ,  b  .  On  noting  that  E  (In  x.)  -  E  (x  )  by  strict  mono- 
n  n  n  l  n  l 

tonicity  of  In,  it  is  apparent  that  the  results  are  the  same  whether  we 

perform  the  (nonlinear)  procedure  for  X  based  on  X  e  on 

In  X  based  on  In  X  e  /5"(A).  The  same  is  true  for  the  linear  procedure 

since  E  (x. )  •  exp(-(x. /a' )”“)  is  transformed  into 
n  i  in 

-ln(-ln(E  (x.)))  -  a  In  x.  -  o  In  s'. 

Hi  1  n 

If  one  were  unsure  as  to  which  (or  any)  domain  of  attraction  is  appropri¬ 
ate,  regressions  corresponding  to  all  the  candidate  extreme  value  distri¬ 
butions  could  be  performed  and  the  adequacy  of  the  fits  evaluated.  For 
instance,  we  could  perform  one  regression  using  y^  ^’s  and  another  using 

in 

Consider  the  specialization  to  $ (A)  and  assume  the  exponential  tall 


assumption  (2.11)  holds.  We  can  use  either  the  basic  linear  or  the 
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nonlinear  regression  procedure  Co  obcain  estimates  a  ,  b  if  n  is  large 

n  n 

enough  for  the  approximation  (2.22)  to  be  good.  Then  without  further 
simulation  or  regression,  we  can  use  (2.11a)  and  (2.11b)  to  obtain 

I 

^  ^  M  ^  A 

estimates  a  ,  -  a  ,  b  ,  ~  a  ln(n'/n)  +  b  .  The  linear  or  nonlinear 
n  n  n  n  n 

regression  procedure  generally  Increases  in  accuracy  with  k/n,  but  the 
relation  (2.22)  deteriorates  as  n  decreases.  Thus  for  a  fixed  amount  of 
simulation  k,  n  should  be  chosen  as  small  as  the  adequacy  of  (2.22) 
allows . 

The  relations  (2.11a)  and  (2.11b)  suggest  a  possible  variation  to 
either  the  basic  linear  or  nonlinear  regression  procedure:  simulate  k 
cycles  and  substitute  aQ  -  a  \  bQ  •  a  1ln(nb),  in  the  rlght~h4nd  side  of 
the  regression.  Now  form  the  empirical  d.f.  as  in  the  basic  procedure  for 
a  fixed  n  -  n^ .  Do  this  for  several  values  of  n ^ ,  pool  the  sample ,  and 
perform  one  regression  to  estimate  a  \  b.  Assuming  (2.11a)  and  (2.11b) 
hold  adequately,  this  variation  extracts  more  information  from  a  given 
amount  of  simulation.  This  comes  at  the  cost  of  increased  complexity,  and 
has  not  been  attempted. 

3.2  Modifications  for  the  Discrete  State  Space  Case 

Now  assume  the  regenerative  process  of  interest  has  a  discrete  state 
space  and  h|  If  the  basic  (continuous  state  spsce)  procedures  were 

used,  there  could  now  be  ties  among  the  y  ’s.  So  we  simulate  k  cycles 

1  f  K 

of  the  process  and  let  ^  k  >  72  k  >•••>  k  be  the  di8tinct  Individual 

cycle  maxima,  occurring  respectively  in  N^  ^.....N^  different  cycles. 

The  number  of  N  's  could  be  quite  small  even  for  k  lasge.  Letting 
1,K  i-1 

N  -  0  and  m  ,  ■  Z  N.  ,  y  is  the  maximum  value  in 

U  p  iC  x  )  K.  .  a  X  ,  K  X  g  1C 


I 


k-m  -jl 
i,k 

n-1 


of  the  (j|)  subsets  of  n  cycles,  for  i  <  JL.  Thus  En(y^  ■  1,  and  if 
for  instance  -  2,  then  EJy^)  -  ^ 

etc.  So  E^(y^  can  still  be  obtained  by  a  simple  recursion.  In 

practice,  we  would  only  consider  yi  ^  for  1  <_  1q,  where 

1^  ■  max{i  :  EnC>ri  >_  *}  for  some  small  «  .  Now  1^  Is  a  function  of 

che  sample  path  since  the  number  of  duplications  is  not  known  a  priori. 

'Sail  the  procedure  where  we  throw  out  y^  ^  and  do  a  regression  as  in  the 

basic  (continuous  state  space  case),  the  basic  full  delete  procedure 

(either  linear  or  nonlinear).  The  basic  partial  delete  procedure  may  not 

completely  throw  out  y^  Instead,  we  consider  a  modified  empirical  d.f. 

as  follows:  if  2  and  say  for  example  *  3,  then  take 

k-n  k-n-1  ,  _  .  .  k-n  k-n-1  .  k-n-2 

- k =r  and  En(y2,k>  “  — - 


etc. 


k  k-1  “““  “nv/2,k'  "  k  k-I  k-2  » 

We  only  alter  E  (y.  .),  leaving  E  (y  ,  ),  i  >  1  as  before.  If 
n  ^  i  k  n  l  >  k. 

N  ■  1,  the  basic  full  and  partial  delete  procedures  are  equivalent; 

1  >  K. 

otherwise  they  differ  only  in  the  use  of  an  additional  regression  data 
point  by  the  partial  delete  procedure.  Most  of  the  variations  mentioned  in 
the  continuous  state  space  case  are  also  available  in  the  discrete  case. 

A  A 

After  having  obtained  estimates  a^,  b^,  we  would  approximate 

f  i  .  A  a  ^  A 

P{  max  HT  <  x;  by  exp(-exp(-a  (x  +  6-b  )))  where  6  represents  a 

l<j<n  3  -  n  n 

continuity  correction  such  as  6  -  0  or  1/2.  Estimating  the  d.f.  of 


max  M.  is  tougher  in  the  discrete  than  the  continuous  case  for  the 
l<j<n  3 


17 


following  reasons:  (i)  a  continuity  correction  (6)  is  needed,  (ii)  true 
convergence  to  the  extreme  value  distribution  is  not  attained,  suggesting 

that  the  approximation  on  which  the  regression  equation  is  base!  may  not  be 

» 

as  good  as  in  the  continuous  case,  and  (iii)  due  to  duplications  associated 

with  the  y  's,  the  number  of  data  points  for  the  regression  is  greatly 

l,k 

reduced  relative  to  the  continuous  case,  for  a  given  k,n,  and  cutoff  e. 

The  basic  procedures  are  illustrated  in  Section  6  for  the  M/M/1 
queue.  Before  leaving  this  section  we  point  out  some  other  relevant 
references.  The  reliability  theory  literature  contains  many  references  on 
the  problem  of  testing  whether  observations  come  from  an  exponential  or 
extreme  value  d.f.  and  of  estimating  the  associated  parameters*  Two  useful 


places  to  find  such  papers  are  EPSTEIN  (1960)  and  MANN,  SCHAFER,  and 
SINGPURWALLA  (1974),  Chapter  5.  PICKANDS  (1975)  has  developed  a  method  for 
determining  which  G  d.f.  is  appropriate  for  a  given  set  of  observations. 
His  method  uses  a  random  and  Increasing  number  of  the  y  's  as  n 

j.n 

increases.  The  method  is  expensive  computationally  and  emphasizes  an 
aspect  of  the  extreme  value  problem  which  is  not  of  great  concern  for 
simulation.  Finally,  WEISSMAN  (1978)  contains  another  method  for  estimat¬ 


ing  the  constants  a  and  b  . 

n  n 


k.  The  GI/G/1  Queue 

The  Gl/G/1  queue  and  the  birth-death  processes  treated  in  Section  5 

are  among  the  very  few  regenerative  processes  for  which  we  know  the 

values  of  a  and  b  .  For  this  reason  these  processes  are  excellent 
n  n 

candidates  for  testing  the  effectiveness  of  the  estimation  procedures 
proposed  in  Section  3* 
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In  the  Gl/G/l  queue  we  assume  customer  0  arrives  at  tQ  =  0, 

finds  .a  free  server,  and  experiences  a  service  time  v^.  Customer  n 

arrives  at  time  t  and  experiences  a  service  time  v  .  Customers  are 
n  n 

served  in  their  order  of  arrival  and  the  server  is  never  idle  if  customers 
are  waiting.  Let  the  interarrival  times  t  -  cn_^  =  un>  n  >  1.  We 
assume  the  two  sequences  (vn:n  ^  0)  and  (un:n  ^  1}  each  consist  of 
i.i.d.  r.v. 's  and  are  themselves  independent.  Let  Efu^)  =  and 
E{ v  ]  =  u*1.  where  0  <  X.  u.  <  ®.  The  traffic  intensity  p  =  Vu  is 

assumed  to  be  less  than  one.  We  exclude  the  deterministic  system  in 

which  both  the  v  's  and  u  '  s  are  degenerate.  Let  the  waiting  time  of 
n  n 

the  nth  customer  be  W  ,  the  workload  ( or  virtual  waiting  time)  at 

n  » 

time  t  be  V  and  the  number  of  customers  in  the  system  at  time  t 

be  Q£.  Also  set  Wr  a  max(W^ :0  <  j  <  n),  V£  =  sup{Vg:0  <  s  <  t), 

* 

and  Q  =  sup{Q  :0  <  s  <  t).  Let  X  =*v  ,  -  u  n  >  1,  and  set 

t  s  “ ’  —  n  n- 1  n'  — 

S  »  X,  +  ♦••  +  X  for  n  >  1  and  S-  »  0.  If  n.  denotes  the  number 
n  1  n  —  0  j 

of  customers  served  in  the  jth  busy  period,  then  n^  is  related  to  the 
partial  sum  process  (S^zn  ^  0}  since 

n.  =  inf(n  >  0:S  <  0)  . 

1  n  — 

When  p  <  1.  we  have  ra  =  E(n,)  <  w.  Also  -S  is  the  length  of  the 
first  idle  period.  We  assume  that  has  an  aperiodic  d.f.  (support 

is  not  concentrated  on  a  set  of  points  of  the  form  0,  +h,  +2h,  4-Jh,  ...), 
that  there  exists  a  positive  number  k  such  that  E(exp(xX^))  *  1, 
and  0  <  =  E(X^  exp(/OC^)j  <  ».  These  assumptions  will  normally  be 

satisfied  if  the  d.f.  of  vQ  has  an  exponentially  decaying  tail;  e.g.. 


19 


r 

i 


when  vQ  has  a  gamma  distribution.  Under  these  conditions  we  know 
(see  IGLEHART  (I972))  that 


(*.l)  . 

a% 

1 

* 

1 

*-• 

In  Cln)/(c"1  =^A^”(x)  . 

and 

(fc.2) 

.  *  -1 
(Vt  -  * 

In  c2t)/x_1  =»  A^m(x)  , 

where 

*Sn,  2 


*v0 

o,  =  E(e  °)Cl  . 


Thus  to  use  (4.1)  and  (4.2)  for  estimating  the  d.f. 's  of  W  and  V 

kS 

we  need  only  estimate  m  and  E{e  *-),  assuming  that  k,  and 

kv  o 

E(e  )  can  be  calculated  numerically.  In  the  special  case  of  M/G/l 

- 1  kS 

queues  no  estimation  is  required,  since  m  =  (l-o)  and  E{e  nl} 

=  \  (A+-k)  .  If  the  simulation  is  carried  out  for  a  fixed  number  of  cycles, 
then  counterparts  of  (4.1)  and  (4.2)  hold  with  the  exponents  of  A 
removed. 

The  queue-length  process  {Qc : t  >  0}  is  discrete-valued  and  the 
associated  d.f.  of  M*  €  <2.  Hence  a  limit  theorem  comparable  to  (4.1) 
or  (4.2)  does  not  exist.  Instead  we  must  seek  results  like  (2.12)  and 
(2.13)  °r  (2.15)  and  (2.16).  Unfortunately,  these  results  are  only 
known  fcr  the  M/G/l  and  GI/M/1  queues;  see  COHEN  (1969),  Theorems  7.2 
and  7.5.  Let  M.*  =  supCQ^iTj  ^  <  t  <  Tj ) .  Then  for  an  M/G/l  queue 
the  counterpart  of  (2.  16)  is 
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(*•3) 

where 


P{  max  M*  <  1) 
1  <  j  <  n  J 


-a(  *-bn) 

exp(-e 


) 


a  =  In  ( {  'K+k)  /  ~K) 

and 

b  =  a"  ^  In  (  c0n) 
n  * 

On  the  other  hand,  for  GI/M/ 1  queues  (4.3)  holds  with  a  =  ln( (u-k) /u) 

and  the  same  value  for  b  .  Tables  1  and  2  contain  the  values  of  m.  k. 

n  *  * 

c^,  and  Cg  for  the  M/M/1  and  M/Eg/ l  queues  as  a  function  of 
the  traffic  intensity  p. 


(4.4)  EXAMPLE.  M/M/1  queue .  For  this  queue  we  are  able  to  calculate 
exactly  the  distribution  of  associated  with  the  waiting  time  process, 

the  virtual  waiting  time  process,  and  the  queue  length  process.  The  tall 
of  the  distribution  of  M*  for  the  waiting  time  process  is  given  in  COHEN 
(1969),  p.  606,  or  can  be  calculated  directly  from  results  in  IGLEHART 
(1972),  p.  630.  The  result  is 

-kx 

P{m[  >  x}  -  P(1~2''-kx-  ~  P(l-P)e~kx  as  x  +  -  . 

1-p  e 

The  corresponding  ''suit  for  the  virtual  waiting  time  process  is  given  in 
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TABLE  1 


Parameter  Values  for  M/M/l  Queue  with  4  =  10 


p 

m 

K 

°i 

c2 

1 

1.11 

9 

.900 

.09 

•  9 

2 

1.25 

8 

.400 

.  16 

.8 

3 

1.43 

7 

.233 

.21 

•7 

4 

1.67 

6 

.150 

.21+ 

.6 

5 

2.00 

5 

.  100 

.25 

•  5 

6 

2.50 

4 

.067 

.21+ 

.4 

7 

3.33 

3 

.01+3 

.21 

•  3 

3 

5.00 

2 

.025 

.16 

.2 

9 

10.00 

1 

.011 

.09 

.1 

95 

20.00 

0.5 

.005 

.0475 

.05 

99 

100.00 

0.  1 

.001 

.0099 

.01 

TABLE  2 


Parameter  Values 

for 

M/Eg/l 

Queue  with 

H  =  10 

P 

m 

K 

•■v 

C1 

c2 

.1 

i.ii 

15.00 

.3375 

.1562 

2.5 

.2 

1.25 

12.60 

.2016 

.231+6 

1.7119 

.3 

1.43 

10.61 

.1395 

.2874 

1.3038 

.1+ 

1.67 

8.83 

.1012 

.3179 

1.0201 

.5 

2.00 

7.19 

.0741 

.3263 

0-7957 

.6 

2.50 

5.64 

.0534 

.3118 

0.6050 

.7 

3.33 

4.  16 

.0367 

.2733 

0.1+357 

.8 

5.00 

2.73 

.0227 

.2091+ 

0.2809 

•  9 

10.00 

1.35 

.0106 

.1188 

0.1366 

•  95 

20.00 

O.67 

.0051 

.0630 

0.0671+ 

•  99 

100.00 

0.13 

.0010 

.0132 

0.0131+ 
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[4],  p.  606,  or  can  again  be  calculated  from  [10],  p.  632.  The  result  is 


-kx 

PK  >  X>  ’  (1~P)!w  ~  (l“P)e_kX  as  x  ♦  »  . 

1-pe 

The  corresponding  result  for  the  queue  length  process  is  given  in  Example 
(5.4).  ^ 


5 .  Birth-Death  Processes 

A  second  class  of  regenerative  processes  for  which  theoretical 
results  are  available  is  birth-death  processes  in  discrete  or  continuous 
time.  Let  (X^in  >  0}  be  a  discrete  time  Markov  chain  with  state  space 
E  =  (0,  1,  2,  ...]  and  transition  probabilities  given  by 

IV  J  *  ‘-1 

(5.1)  Pt)  -  <  Pt,  j  »  1+1 

'  0  ,  other  j, 

where  =  0,  PQ  =  ^  and  the  other  q^'s  ant^  p^'s  are  positive. 

This  chain  will  automatically  be  both  irreducible  and  periodic.  Further¬ 
more,  recall  that  it  will  be  recurrent  if  and  only  if 

00 

Z  (^T.  Px)  1  =  » 
j=i  i  J 

where  =  1  and  n .  =  ( p_  • • •  p.  •••  q,).  We  assume  the  chain 

0  J  0 

is  recurrent.  It  will  be  positive  recurrent  if  and  only  if 
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] 


00 


Z  it  .  <  oo  . 

j=0  J 

Next  define 


T^(k)  =  inf {n  >  OiX^  =  k),  k  E  E  , 


the  first  entrance  time  to  state  k.  Let  [ • )  =  P{*!XQ  =  i},  the 
conditional  probability  of  an  event,  given  XQ  =  i.  Then  our  concern 
here  will  be  in  the  probability,  given  X^  =  i,  of  the  Markov  chain 
entering  state  n  before  it  enters  state  0.  Let  this  probability  be 
denoted  by 


rjn)  =  P.(T1(n)  <  t]_(0))  ,  i  €  { 1,  2,  ...,  n-1)  . 


Fortunately,  this  probability  has  been  calculated  and  in  particular 

r0(°)  =  r1(n)  =  (1  +  S  (rc.p.)'1)'1  > 

i=  1 

see  CHUNG  (I960),  p.  68.  Note  that  lim  r^n)  =  0  when  the  chain  is 

n  — » 00 

recurrent,  in  keeping  with  our  intuition.  Define 

=  sup{Xn:0  <  n  <  t^(0)  -  1}  . 

Then 

+  /  n  .iv- 1 

(5.2)  P0(M1  >  n]  =  rQ(n^l)  =  I  Z  ( it  p  )  )  ,  as  n  -»»  . 

'  i=0  / 
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Suppose  now  that  we  have  a  birth-death  process  (X t : t  >  0}:  a 
continuous  time  Markov  chain  with  state  space  E  =  {0,  1,  . ..)  and 
embedded  jump  chain  whose  probabilities  are  given  by  (5.1).  As  above, 

I 

define  the  first  entrance  time  to  state  k  and  the  maximum  in  the  first 
cycle  by 


TL(k)  =  inf  ( s  >  0:Xs_  ±  k,  Xg  =  k) 

and 

M*  =  sup(Xt:  0  <  t  <  t^O)}  . 

Because  of  the  path  structure  of  the  birth-death  process,  it  is  easy  to 
see  from  (5.2)  that 

(5-3)  VM1  >  n]  =  [l  +  \  (*i  V"1]  ' 

where  V  [resp.  are  the  birth  [resp.  death]  parameters  and 

jIq  =  1,  ni  =  ^i)  *  The  same  argun,ant  can  of 

course  be  used  to  show  that  (5-3)  also  holds  for  semi-Markov  processes 
with  embedded  jump' chain  whose  probabilities  are  given  by  (5-1). 


(5.L)  EXAMPLE.  M/M/s  queue.  The  queue- length  process, 
is  a  birth-death  process  with  parameters  =  X  and  ^  = 
j  >  0.  Assume  the  queue  has  traffic  intensity  p  =  "h/\xs  < 
and  sufficient  condition  for  recurrence.  Then  from  (5.3) 


{Qt:t  >  0), 
n  (j  A  s), 

1,  a  necessary 
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VMt 

>  n}  =  | 

'  i  .-/r1 

■i=0  1  J 

-1 

■  s  .  e,  n-( s+1) 

s  <x>  i:  ♦tHt!  1 

■i=0  *  s8  ps+1  i=0 

Asymptotically^ 

as  n  - 

-4  00 

(  (s3/s.')n_1  ,  P  = 

foK 

>  n)  ~ 

) 

l  (sS(l-p)/s.')pn  ,  p< 

Thus  for  p  <  1 

we  can 

use  (2.16)  to  obtain 

P  {  max  Mt  <  ]  =  exp( -e  bn^)  , 

0  1  <  j  <  n  2 


where  a  =  In  p  and  =  a  In  (n  s  (l-p)/s,).  Note  that  Phis  is 
consistent  with  ( i+ . 3)  -  ^ 
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6.  Numerical  Results 


A  simulation  of  the  M/M/1  queue  was  performed  to  help  assess  the 

effectiveness  of  the  basic  procedures  proposed  in  Section  3.  Our  goal  is 

to  estimate  the  d.f.  of  max  M+,  where  M*  is  either  the  maximum  waiting 

l<j<n  S  S 

time  W,  virtual  waiting  V,  or  number  of  customers  Q  in  the  system  in 

cycle  j.  The  processes  {W^  :  n  ^  0}  and  {V  :  t  ^  0}  have  continuous 

state  spaces,  {Qt  :  t  ^  0}  has  a  discrete  state  space,  and  A  is  the 

appropriate  extreme  value  distribution  in  all  three  cases.  Using  estimates 

a  ,  b  ,  we  approximate  P{  max  M*  <  x}  by  A((x  +  6  -b  )/a  )  where  6*0 
n  n  Kj<n  J  —  n  n 

for  and  V  ,  and  we  try  6-0  and  6  -  1/2  for  Qt»  Theoretical 

values  of  a  and  b  as  well  as  the  exact  values  of  P{  max  -  M*  <  x} 
n  n  K  j<n  i  ~ 

are  available  from  the  results  in  Sections  4  and  5.  All  notation  and 

conventions  in  this  section  are  as  in  Section  3. 

The  random  number  generator  used  was  the  DEC-20  FORTRAN  "RAN" 

function.  In  the  terminology  of  Section  3,  all  experimental  results 

reported  are  for  basic  linear  or  nonlinear  regression  procedures.  Several 

values  of  k,n,  and  traffic  intensity  p  were  tested.  Both  the  linear  and 

nonlinear  regression  procedures  were  tried  on  Wq,  while  only  the  linear 

regression  procedure  was  tried  on  V  and  Q^.  The  nonlinear  regression 

was  performed  using  SUBROUTINE  LSQFDN  of  the  National  Physical  Laboratory 

Algorithms  Library.  Both  full  and  partial  delete  as  well  as  continuity 

corrections  of  6  -  0  and  6  -  1/2  were  tried  on  Q^.  For  Qt»  a  cutoff 

of  £  —  .0001  was  used.  For  W  and  V  ,  a  cutoff  of  £  -  .001  was 

n  t 

used.  Thus  for  W  and  V  .  i.  is  as  follows: 
n  t  0 


28 


k 


n  *0 


100,000  2000  342 

100,000  1000  686 

50,000  1000  341 

20,000  1000  135 

50,000  250  801* 

25,000  250  678 

10,000  250  270 


The  starred  entry  Indicates  that  for  k  «  50,000,  n  ■  250,  i^  corresponds 

to  e  «  .01755  rather  than  to  e  m  .001.  Note  that  in  all  cases,  there 

were  sufficient  points  for  the  regression  for  Q  ;  and  y  ,...,y 

t  1 f K  10 

were  always  distinct  for  and  V  . 

Tables  3-5  contain  the  results  of  the  simulation  for  estimating  a 

n 

and  b  .  Tables  6-10  contain  the  estimates  of  P{  max  M*  <  x}  by  • 
n  1<  j<n  i  ~ 

A((x  +  6  -  b  )/a  ),  using  the  estimated  values  of  a  ,  b  shown  in  Tables 
n  n  n  n 

3-5.  The  entries  contained  in  the  tables  are  the  sample  means  of  the 
various  estimates  over  the  number  of  replications  and  the  half-length  of  a 
symmetric  90%  confidence  interval  about  the  sample  mean.  For  example,  in 
Table  4,  take  k  =  50,000,  n  ■  1000,  and  50  replications.  Then  a  90%  con¬ 
fidence  interval  for  a100Q  based  on  50  replications  of  the  linear 
regression  procedure  is  [.1966  -  .0039,  .1966  +  .0039].  The  corresponding 
true  value  of  a^Q00  *  *^*  Tables  6-10  report  the  true  values  of 

A((x  -  b  )/a  )  and  P{  max  M*  <  x}  for  various  x,  as  well  as  the 
«  n  l<j<n  J  - 

corresponding  values  of  A((x  +  6  -  b  )/a  )  using  estimated  a  ,  b  .  For 

an  n  n 

example,  in  Table  9,  take  k  ■  50,000,  n  *  250,  and  100  replications.  Then 

the  true  values  of  A((9  -  b,e-)/a_cn)  and  P{  max  M*  <  9} 

250'  250'  l<j<250  J  ~ 

are  respectively  .7834  and  .7831,  while  the  sample  mean  of 

A((9.0  -  b  )/a  )  is  .7832  with  a  corresponding  90%  confidence  interval 
n  n 

half-length  of  .0062. 
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TABLE  3 


Estimates  of  an  and  bn  for  {*n  :  n  >_0}  In  the  M/M/1  queue  with  |1  ■  10 


k/n/#repl . 


100,000/1000/50  .5 


50,000/1000/50  -5 


20,000/1000/50  .3 


50,000/250/100  .5 


25,000/250/100  .5 


10,000/250/100  .5 


100,000/1000/50 


50,000/1000/50 


20,000/1000/50 


100,000/2000/50 


true  values  linear  regression  I  nonlinear  regression 


.2  1.1043 


.2  1.1043 


.2  1.1043 


.2  .8270 


.2  .8270 


.2  .8270 


4.4998 


4.4998 


4.4998 


5.1930 


TABLE  4 


Estimates  of  an  and  bn  for  {Vf  :  t  >_ 0}  In  the  M/M/1  queue  vith  |i  ■  10 


true  values 

linear  regression 

k/n//rep I . 

P 

*n 

*n 

* 

•n 

* 

100,000/1000/50 

.5 

.2 

1.2429 

.1980 

1.2451 

.0021 

.0034 

50,000/1000/50 

.5 

.2 

1.2429 

.  t966 

1.2422 

.0039 

.0057 

20,000/1000/50 

.5 

.2 

1.2429 

.1917 

1.2371 

.0062 

.0093 

50,000/250/100 

.5 

.2 

.9657 

.1977 

.9692 

.0016 

.0018 

25,000/250/100 

.5 

.2 

.9657 

.1966 

.9673 

.0018 

.0029 

10,000/250/100 

.5 

.2 

.9657 

.1946 

.9648 

.0027 

.0047 

100.000/1000/50 

.9 

1.0 

4.6052 

.9713 

4.6050 

i 

I 

.0115 

.0214 

50,000/1000/50 

.9 

1.0 

4.6052 

|  .9819 

4.6308 

.0163 

.0276 

20,000/1000/50 

!  •» 

1.0 

|  4.6052 

.9758 

4.6401 

.0251 

.0442 

100,000/2000/50 

i 

•9 

1.0 

5.0695 

.9763 

5.2840 

.0161 

.0293 
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TABLE  5 


Estimates  of  an  and  bn  for  {Q+  :  t  >_  0}  In  the  M/M/I  queue  with  |i  »  10 


« 

full  delete 

partial 

1  delete 

true  values 

1  linear  regression 

linear  regression 

k/n/#rep 1 . 

P 

•n 

*>n 

an 

»n 

a* 

«n 

*e 

«>n 

100,000/1000/50 

.5 

1 .4427 

8.9658 

1.4350 

8.9779 

|Pl 

8.9763 

.0409 

.0394 

.0395 

50,000/1000/50 

.5 

1.4427 

8.9658 

1.4224 

8.9438 

8.9426 

.0519 

.0536 

.0540 

20,000/1000/50 

.5 

1.4427 

8.9658 

1.3783 

8.8994 

1.3801 

8.8980 

.0611 

.0771 

.0609 

.0758 

50,000/250/100 

.5 

1.4427 

6.9658 

1.4347 

6.9853 

1.4336 

6.9842 

.0253 

.0179 

.0212 

'  .0180 

25,000/250/100 

.5 

1.4427 

6.9658 

1.4224 

6.9725 

1.4179 

6.9695 

.0288 

.0252 

.0301 

.0257 

10,000/250/100 

.5 

1 .4427 

6.9658 

1 .4039 

6.9491 

1.3986 

6.9437 

i 

.0345 

.0379 

.0348 

.0374 

100,000/1000/50 

.9 

9.4912 

43.7087 

9.2227 

43.7424 

9.2250 

43.7421 

.1945 

.2232 

.1921 

.2226 

50,000/1000/50 

.9 

9.4912 

43.7087 

i 

9.3229 

43.9337 

9.3203 

43.9316 

.2351 

.3107 

.2365 

.3116 

20,000/1000/50 

.9 

9.4912 

43.7087 

9.1984 

43.8827 

9.2080 

43.8893 

.3155 

.4561 

.3178 

.4580 

100,000/2000/50 

.9 

9.4912 

50.2875 

9.2249 

50.1689 

9.2261 

50.1684 

i 

.2333 

.3108 

.2297 

.3085 
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TABLE  6 


Estimate*  of  P(mk  H*  <  x>  tor  {#  :  n  >  0}  In  the  M/M/1  queue  with  4  ■  10,  p  ■  .5 

l<j  <n  J  ~  n  — 


regression 

method 

h/n/#r*pl. 

A<(x-bn)/an)/x 

P(  MX  M*  <  x) 

1 <n  J  — 

.25/1 .0390 
.2493 

.50/1.1776  .75/1.3535 

.4996  .7499 

.90/1 .5544 
.9000 

.99/2.0243 

.9900 

1 1  near 

100,000/1000/50 

.2457 

.4971 

.7489 

.8997 

.9900 

.0052 

.0072 

.0062 

.0037 

.0007 

non  1 i near 

100.000/1000/50 

.2442 

.4966 

.7492 

.8999 

.9899 

.0054 

.0079 

.0074 

.0045 

.0008 

1 1  near 

50,01/0/1000/50 

.2511 

.5040 

.7535 

.9015 

.9899 

.0082 

.0117 

.0104 

.0062 

.0012 

non  1 1  near 

50,000/1000/50 

.2489 

.5035 

.7539 

.9014 

.9896 

.0071 

.0121 

.0116 

.0074 

.0015 

1 1  nea" 

20,000/1000/50 

.2626 

.5193 

.7638 

.9053 

• 

.9897 

.0136 

.0194 

.0174 

.0107 

.0022 

non  1 1  near 

20,000/1000/50 

.2576 

.5208 

.7683 

.9079 

.9898 

.0131 

.0188 

.0174 

.0110 

.0023 

A((x  -b„)/an)/x 

P(  max  M*  <  x) 

1<J  <n  J  — 

.25/. 7617 

.50/. 9003  .75/1.0762 

.90/1.2771 

.99/1.7471 

.2471 

.4986 

.7496 

.8999 

.9900 

1 1  near 

50,000/250/100 

|  .2438 

.4963 

.7493 

.9003 

.9901 

.0027 

.0042 

.0039 

.0023 

.0005 

non  1 1  near 

50.000/250/100 

.2432 

.4951 

.7481 

.8996 

.9900 

.002B 

.0040 

.0038 

.0023 

.0005 

1 1  near 

25,000/250/100 

.2464 

.5002 

.7525 

.9020 

.9903 

.0042 

.0060 

.0052 

.0030 

.0005 

non  1 1  near 

25,000/250/100 

.2454 

.5001 

.7531 

.9025 

.9904 

.0045 

.0060 

.0053 

.0032 

.0005 

1 1  near 

10,000/250/100 

.2529 

.5092 

.7593 

.9052 

.9906 

.0070 

.0097 

.0082 

.0049 

.0008 

non  1 1  near 

10,000/250/100 

.2494 

.5087 

.7615 

.9071 

.9909 

.0076 

.0098 

.0081 

.0043 

.0008 

33 


TABLE  7 


Estimate*  of  P(max  M*  <  x)  for  {w  :  n  >  0}  In  the  M/M/1  queue  »lth  u  «  10,  p  ■  .9 
1  <n  J  —  «  — 


regression 

method 

k/n/#rep 1 • 

At(x-bn)/an)/x 

P(  max  M|  <  x) 

1  <J  <n  J  — 

.25/4.1732 

.2454 

.50/4.8663  .75/5.7457 

.4977  .7494 

.90/6.7502 

.8999 

.99/9.1000 

.9900 

1 1  near 

100,000/1000/50 

.2484 

.5049  .7577 

.9054 

.9910 

.0065 

.0089  .0074 

.0042 

.0007 

non  1 1  near 

100,000/1000/50 

.2492 

.5050  .7569 

.9045 

.9907 

.0065 

.0089  .0079 

.0049 

.0008 

1  i  near 

50,000/1000/50 

.2415 

.4952  .7490 

.9001 

.9900 

.0079 

.0116  .0101 

.0059 

.0010 

non  1  inear 

50,000/1000/50 

.2446 

.4990  .7514 

.9006 

.9897 

.0082 

.0117  .0111 

.0069 

.0012 

1  i  near 

20,000/1000/50 

.2422 

.4968  .7491 

.8991 

.9894 

.0131 

.0184  .0161 

.0094 

.0017 

non  1 1  near 

20,000/1000/50 

.2442 

.5045  .7558 

.9012 

.9888 

.0132 

.0201  .0191 

.0119 

.0023 

A((x-bn)/an) 

P<  max  M'*'  <  x> 

1<J <n  J  ~ 

.25/4.8663 

.50/5.5595  .75/6.4389 

.90/7.4433 

.99/9.7931 

.2477 

.4989  .7497 

.9000 

.9900 

1 1  near 

100,000/2000/50 

•  .2552 

.5100  .7590 

.9049 

.9906 

.0090 

.0121  .0101 

.0059 

.0010 

non  1 1  near 

100,000/2000/50 

.2539 

.5112  .7607 

.9055 

.9905 

.0092 

.0132  .0119 

.0070 

.0012 
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TABLE  8 

Estimates  of  Ptmax  M*  <  x)  for  {v  :  t  >  0}  In  the  M/M/1  queue  xlth  |i  •  10, 
1  <J  <n  J  ~  *  ~ 

using  linear  regression 


A((x-bn)/anl/x 

P<  max  M*  <  x) 

1<J  <n  J  ~ 

k/n/#repl. 

P 

.25/1.1776 

.2493 

.50/1.3162  .75/1.4921 

.4996  .7499 

.90/1 .6930 
.9000 

.99/2.1630 

.9900 

100,000/1000/50 

.9006 

.9901 

I 

.0037 

.0007 

50,000/1000/50 

I 

.2523 

.5062 

.7555 

.9025 

.9900 

■  I 

.0082 

.0122 

.0111 

.0067 

.0012 

20,000/1000/50 

| 

.2631 

.5224 

.7676 

.9077 

.9902 

H 

.0134 

.0196 

.0174 

.0106 

.0020 

At  Cx-bn./an)/x 

Pt  max  M*  <  x) 

1  <J  <n  J  ~ 

.25/. 9003 

.50/1.0390  .75/1.2148 

.90/1.4157 

.99/1.8857 

.2471 

.4986 

.7496 

.8999 

.9900 

50,000/250/100 

n 

.2430 

.4959 

.7494 

.9005 

.9901 

1 1 

.0027 

.0040 

.0039 

.0025 

.0005 

25,000/250/100 

| 

.2466 

.5005 

.7528 

.9021 

.9904 

| 

.0042 

.0062 

.0054 

.0032 

.0005 

10.000/250/100 

| 

.2518 

.5076 

.7579 

.9044 

.9905 

Bl 

.0070 

.0097 

.0082 

.0049 

.0008 

At  tx-bn)/an)/x 

P(  max  M|  <  x) 

1<J  <n  J 

.25/4.2785 

.50/4.9717  .75/5.8511 

.90/6.8555 

.99/9.2053 

.2454 

.4977 

.7494 

.8999 

.9900 

100,000/1000/50 

mm 

.2481 

.5048 

.7576 

.9054 

.9910 

| 

.0065 

.0089 

.0074 

.0042 

.0007 

50,000/1000/50 

| 

.2414 

.4953 

.7491 

.9002 

.9900 

| 

.0079 

.0116 

.0100 

.0059 

.0010 

20,000/1000/50 

D 

.2414 

.4958 

.7484 

.8987 

.9893 

H 

.0129 

.0184 

.0159 

.0094 

.0017 

At<x-bn)/an)/x 

P(  max  M*  x) 

1  <J <n  J 

.25/4.9717 

.50/5.6648  .75/6.5442 

.90/7.5487 

.99/9.8985 

.2477 

.4989 

.7497 

.9000 

.9900 

100,000/2000/50 

.9 

.2550 

.5100 

.7591 

.9050 

.9906 

.0090 

.0117 

.0102 

.0059 

.0010 
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TABLE  9 


Estimates  of  Pfmax  M*  _<_x>  for  {Q  :  t  0}  In  the  M/M/1  queue  with  p  ■  10, 
1<J<n  J  f 

p  »  .9,  using  full  delete  linear  regression,  with  continuity  correction  6  ■  0 


A(  {x-bn)/an 

>/x 

P(  max  M*  < 
1  <J  <n  J 

x) 

k/n//repl . 

.3766/9 

.3761 

.6137/10 

.6135 

.7834/11 

.7833 

.9408/13 

.9408 

.9924/16 

.9924 

100,000/1000/50 

.3769 

.0106 

.6153 

.0129 

.7839 

.0114 

.9394 

.0057 

.9915 

.0013 

50,000/1000/50 

.3890 

.0144 

.6264 

.0173 

.7908 

.0151 

.9407 

.0077 

.9913 

•0020 

20,000/1000/50 

.4073 

.0213 

.6447 

.0240 

.8032 

.0198 

.9443 

.0092 

•9916 

•0022 

AU*-bn)/an 

)/x 

P(  max  M*  < 
1<J<n  J 

x> 

.3766/7 

.3744 

.6137/8 

.6128 

.7834/9 

.7831 

.9408/1 1 

.9408 

.9924/14 

.9924 

50,000/250/100 

.3732 

.0048 

.6130 

.0065 

.7832 

.0062 

.9399 

•  0032 

.9918 

•0008 

25,000/250/100 

.3784 

.0066 

.6186 

.0088 

.7872 

.0080 

.9410 

.0040 

.9918 

.0010 

10,000/250/100 

.3870 

.0101 

.6273 

.0118 

.7933 

.0101 

.9427 

.0050 

.9919 

.0013 
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TABLE  10 


Estimates  of  P(imx  M  x)  for  {p  :  t  0}  In  the  M/M/1  queue  with  |i  «  10, 
1<j<n  J  + 

p  ■  >9,  using  full  delete  linear  regression,  with  continuity  correction  6  “  0 


k/n/#rep 1 . 

At  (x-bn)/an)/x 

P(  max  M*  <  x) 

1<J  <n  J  ~ 

.2644/41 

.2599 

.5293/48 

.5272 

.7604/56 

.7599 

.9089/66 

.9088 

.9906/88 

.9906 

100,000/1000/50 

.2620 

.5346 

.7678 

.9133 

.9911 

.0070 

.0106 

.0097 

.0059 

.0012 

50,000/1000/50 

.2578 

.5275 

.7604 

.9084 

.9901 

.0095 

.0144 

.0129 

.0074 

.0013 

20,000/1000/50 

.2628 

.5345 

.7647 

.9095 

.9898 

.0147 

.0209 

.0179 

.0104 

.0020 

A((x-b„)/an>/x 

Pt  max  M‘1’  < 

x) 

ljy<n  J 

.2801/48 

.5085/54 

.7695/63 

.9035/72 

.9901/94 

.2779 

.5074 

.7693 

.9034 

.9900 

100,000/2000/50 

.2861 

.5204 

.7796 

.9085 

.9903 

.0104 

.0142 

.0124 

.0077 

.0015 
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Here  are  some  general  observations  on  Tables  3-10.  The  approximation 

of  P{  max  M*  <  x}  by  the  exact  value  of  A((x  -  b  )/a  )  improves  with 
1<  j<n  J  ~  n  n 

greater  x,  greater  n,  and  lesser  p,  and  is  about  equally  good  for  W^, 
Vt#  and  'Qt.  Note  that  confidence  intervals  generally  cover  the  true 
values  (excepting  the  use  of  continuity  correction  6  »  1/2  for  Q^) 
especially  in  the  percentile  evaluations  of  Tables  6-10.  The  bias  and  con¬ 
fidence  interval  widths  tend  to  decrease  as  k  increases,  and  decrease  as 
x  increases  in  the  percentile  evaluations.  The  results  for  p  -  .9  are 
just  about  as  good  as  for  p  *  .5.  The  results  for  are  almost  as  good 

as  for  W  and  V  ,  except  that  the  confidence  interval  widths  are  greater 
n  t 

for  ,  which  can  be  largely  attributed  to  the  greatly  reduced  number  of 
regression  data  points. 

For  the  W  process,  the  linear  regression  procedure  tends  to  produce 
n 

estimates  with  slightly  smaller  biases  and  confidence  intervals  than  does 

the  nonlinear  regression  procedure.  Although  neither  is  a  clear  winner,  in 

this  case  the  linear  regression  would  appear  to  be  better,  if  just  for  the 

reason  of  the  greater  ease  and  lesser  expense  it  entails.  (However,  ;he 

nonlinear  regression  problem  under  consideration  is  a  nice  one,  and  the 

linear  regression  procedure  could  certainly  provide  good  starting  values  of 

a  ,  b  for  the  nonlinear  regression.) 
n  n 

For  the  process,  the  (null)  continuity  correction  of  6-0  is 

superior  to  6  -  .5;  and  thus  results  for  6  -  .5  are  not  shown.  This  is 

no  surprise  in  the  current  case  since  the  approximation  of  tne  exact  d.f. 

of  max  M"t  by  the  value  of  A  using  exact  a  ,  b  has  the  same  error 
l<j<n  J  n’  n 

patterns  for  as  for  and  V£.  The  full  and  partial  delete  pro- 


cedures  give  comparable  results,  so  only  the  estimates  of  a  ,  b  are 

n  n 

shown  for  the  partial  delete  procedure. 

Some  experimentation  on  the  process  (not  reported  in  the  tables) 

using  linear  regression  estimates  a  ,  b  to  arrive  at  estimates  a  .  ■  a  . 

n  n  n'  n 

b^,  -  ln(n'/n)  +  b^  and  corresponding  percentile  evaluations  for  n’ 
cycles  was  performed  using  n'  »  1000  and  n  *  250  and  100.  The  general 
pattern  of  results  is  that  for  fixed  k,  smaller  values  of  n  give  smaller 
confidence  intervals  but  greater  biases,  unless  k  is  very  small.  The  de¬ 
crease  in  confidence  interval  widths  is  mostly  attributable  to  the  use  of 
more  regression  data  points  for  smaller  n.  Again  considering  the  results 
reported  in  the  tables  for  n  ■  1000,  in  many  instances  k  -  50,000 

# 

performs  better  than  k  »  100,000.  All  this  suggests  that  increasing  the 
k/n  ratio  beyond  a  certain  point  may  not  be  worthwhile.  However,  optimal 
values  of  k  and  n  are  very  problem  dependent. 

Most  of  the  cost  in  using  these  procedures  derives  from  the  develop¬ 
ment  and  running  of  the  simulation  model.  Another  substantial  portion  of 

the  cost  is  keeping  track  of  the  y  's.  However  this  can  be  done  once, 

i  >  k 

and  then  empirical  distribution  functions  corresponding  to  the  different 
values  of  n  can  easily  be  constructed.  Different  variations  of  the  re¬ 
gression  procedures  could  be  performed  for  each  value  of  u.  In  short,  the 
variations  are  not  competitors,  but  rather  can  be  used  to  help  assess  the 
validity  of  the  results.  Any  special  insight  to  the  problem  at  hand  could 
also  be  used  to  advantage. 
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